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ABSTRACT 

Context. High-precision cosmology requires the analysis of large-scale surveys in 3D spherical coordinates, i.e. spherical 
Fourier-Bessel decomposition. Current methods are insufficient for future data-sets from wide-field cosmology surveys. 
Aims. The aim of this paper is to present a public code for fast spherical Fourier-Bessel decomposition that can be 
applied to cosmological data or 3D data in spherical coordinates in other scientific fields. 

Methods. We present an equivalent formulation of the spherical Fourier-Bessel decomposition that separates radial and 
tangential calculations. We propose the use of the existing pixelisation scheme HEALPix for a rapid calculation of the 
tangential modes. 

Results. 3DEX (3D Expansions) is a public code for fast spherical Fourier-Bessel decomposition of 3D all-sky surveys 
that takes advantage of HEALPix for the calculation of tangential modes. We perform tests on very large simulations 
and we compare the precision and computation time of our method with an optimised implementation of the spherical 
Fourier-Bessel original formulation. For surveys with millions of galaxies, computation time is reduced by a factor 4-12 
depending on the desired scales and accuracy. The formulation is also suitable for pre-calculations and external storage 
of the spherical harmonics, which allows for additional speed improvements. The 3DEX code can accommodate data 
with masked regions of missing data. 3DEX can also be used in other disciplines, where 3D data are to be analysed in 
spherical coordinates. The code and documentation can be downloaded at http://ixkael.com/blog/3dex 
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1. Introduction 

In the last few decades, cosmology has become a data- 
driven field, where high-precision meas urements of the cos - 
mic microwave bac kground (CMB, e.g.,|Larson et al.|2011[ ) , 
'Schrabback et al7||201Q[ ) and galaxy sur- 



weak lensing 
(e.g. 



Percival et al. 2007 bp have permitted the es- 



veys 

tablishment of a standard cosmological model in which the 
Universe is composed of 4% baryons, 22% dark matter and 
74% dark energy. Some major questions remain, the nature 
of dark matter and dark energy in particular is still not un- 
derstood. Similarly, the initial conditions of the Universe 
are yet to be established and alternative models of gravity 
are still to be tested in comparison with Einstein's general 
relativity. 

New surveys are underway with these science objec- 
tives, e.g. Planck for the CMB (The Planck Colla boration 
[2006]), DES (Dark Energy Survey, |Annis et al.|2005D, BOSS 



[2QC 

(Bar y on Oscillation Spectroscopic Survey, |Schlegel et ah 
2QQ7D, ^ ~ T 
LSS# 

201 Op for "weak lensing and the study of large-scale structure 



gel 

^u. LSST (Large Sy noptic Survey Teles cope, Tyson fc" 
LSS T|2004P and Euclid flLaureijs et al.l2011| [Refregi er et al. 
2010) for weak lensing and the study of large-scale si 
with galaxy surveys. In order to be beneficial, cosmological 
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studies of these surveys need to use high-precision statisti- 
cal methods, such as a full 3D analysis on the sky where 
all- sky 3D surveys are available. 

Several tools have been developed to analyse data on 
the sphere, wh ich is required for a 2D spherical harmonic 
CMB analysis ([Crittenden fc Turok| 19981 |Crittenden|2000 
|G6rski et al.||2002[ |2005| |Doroshkevich et al.||2008p . Weak 
lensing and galaxy survey data can also be analysed tomo- 
graphically (i.e. in 2D slices), but unlike for the CMB, a 
full 3D spherical Fourier-Bessel analysis can also be sought 



dFisher et al . || 199 5 [ [Heavens fc Taylor||1995l [Heave ns 2003; 
|Castro et al.|[2QQ5| jRassat fc Refregier||201lj ). Previous 3D 



data analyses were on relatively small data sets ([Fisher 
[ et al. 119951 |Heaveris & Taylor|1995[[Erdogdu (b) et al.|2006| 
Erdogd u (a) et al.|2006p , but future surveys like Euclid ana 
LSST will provide surveys with billions of galaxies, making 
previous methods for calculating the 3D spectra unfeasibly 
time-consuming. 

In Section |2.1[ we present the theory behind the 3D 
Fourier-Bessel decomposition for infinite and finite contin- 
uous fields as well as the usual method for a discrete sur- 



vey (e.g. galaxy survey). In section 2.2, we present two ad- 
ditional equivalent formulations of the spherical Fourier- 
Bessel decomposition, one of which is central to the 3DEX 
code. In Section [3j we compare the accuracy and calcula- 
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tion time for the usual method used for calculating Fourier- 
Bessel coefficients and methods with the 3DEX code pre- 
sented in this paper. In Section |4j we describe the 3DEX 
library and give examples of how to use it. In Section [5] 
we present our conclusions. We also include an appendix, 
where we discuss the subtleties of the Fourier-Bessel nor- 
malisation. 



2. Theory 

2.1. The spherical Fourier-Bessel decomposition 

In observational cosmology, spherical coordinates (where 
the observer is at the origin) are a natural choice for the 
analysis of cosmological fields. In this system of coordi- 
nates, eigenfunctions of the Laplacian operator are prod- 
ucts of spherical Bessel functions and spherical harmonics, 
i.e. functions j£(kr)Yi m with eigenvalues — fc 2 . For an ho- 
mogeneous three-dimensional field /(r) = /(r, </>) in a 
flat geometry, the spherical Fourier-Bessel decompo sition 
flFisher et aL]1995| |Heavens||2QQ3| |Castro et aL|2005| ) is 



about the underlying matter density field through their spa- 
cial positions. Note that these tracers are subject to various 
distortions and non-linearities, but these are not the pur- 
pose of this work. In this work we only consider linear or 
quasi-linear scales (£ < 50, k < 0.2hMpc -1 ). 

If the only information about the field / is a list of 
coordinates r p = (r p , P , <j> p ) with p = 1, . . . , N (where N is 
the number of galaxies in the latter example), the survey 
may be considered as a superposition of 3D Dirac deltas 
and each coefficient ftrnn can simply be estimated with a 
sum ([Heavens fc Taylor|[l995l |Fisher et aL 1995 Erdogdu 
|(b) et al.||2006| |Abramo et al.||2010p 



N 



/» = ^ 3 )(r-r p ), 
p=i 

JV 

hmn = ^ ke n je(ke n r p )Y e * m (9 p , <j> p ). 



(6) 
(7) 



P =i 



f(r, 



0> 4) = \j\ j dk S fem(k)kj e (kr)Y em (0, </>), (1) 

tm 



This formulation has been used for the analysis of shal- 
low galaxy surveys such as the IRAS 1.2m J survey [~ 6k 
galaxies, |Strauss et al.| 1992 Fisher et al.||199"5] [Heavens 



with the inverse relation 



fc Taylor||1995| ), and th e 2MRS survey"72MASS Redshift 

1. 



fim(k) 



- I d 3 r f(r,0,4>)kj e (kr)Y; m (0, 

7T , 



(2) 



Surv ey, ^ 45k galaxies, |Huchra et al. 201 lj Erdogdu (b) 
et al.||2006| |Erdokdu (a) et al.||2006p . Sincethe time to cal- 
culate equation u\ is proportional to A^n max (f max + l) 2 /2, 
Equation [7] will oecome highly time-consuming when ap- 
plied to larger surveys or when precise decomposition is 
required (large n max and l max )- 



Note that this decompos ition uses the same nota tion as 



Rassat & Refregier (2011) and Castro et al. ( 2005 ), which 2.2. Three equivalent formulations 



is slightly differ ent from the one used in [Lanusse, Rassat 



& Starck ( 2011[ ). The coefficients may be used to calculate 
the 3D power spectrum C(£, fc), defined by 

{hm{k)ft, m ,{k')) = C(l, k)5 D (k - k')5 u >6 mm >, (3) 
a nai've estimator of which is 

C f .(k) = ^^2\f em (k)\ 2 . (4) 

m 

This can be seen as an extension of the usual 2D power 
spectrum (femfe' m '} = Qfe^mm'- The latter arises from 
the spherical harmonic transform of a 2D field given on the 
sphere f(0, (j)) = Y,tm hwYim{0, (j)). 

In practice, surveys will only cover a finite amount 
of volume, limiting the analysis to a sphere of radius R. 
These boundary conditions lead to a discrete spectrum 
{kin}-, which is detailed in the appendices. In this paper, 
we assumed as a boundary condition that / vanishes at 
r = R. The spherical Fourier-Bessel decompo sition becomes 
( |Erdogdu (b) et al.|2006| |Fisher et al]|1995 ) 



In spherical coordinates, since 3D space can be viewed 
as an infinite series of closed shells f2(r), the spherical 
Fourier-Bessel decomposition may also arise from repeated 
2D spherical harmonic t ransforms to which s pherical Bessel 
transforms are applied ( Abramo et al.||2010| ). Formally, the 
field / given on each shell il(r) is first expanded into spher- 
ical harmonics 



f(r,0,4>) = J2fem(r)Y £m (9,ct>), 



(8) 



for which the inversion formula gives harmonics coefficients 
fim(r) depending on the radius r 



hm{r)= / dn f(r,0,4>)Y e * m (e,4>). 
Jn(r) 



(9) 



It is then possible to perform a spherical Bessel transform 



(10) 



fim(r) = \j - / dk fim(k)kji(kr), 



/(r, 6, (fy = ^2 K£nf£rn(k£ n )ke n jt(ke n r)Y £m (6, <j>), (5) 

Imn 

which is exact if the ranges of £,m and n are infinite. The 
Fourier-Bessel coefficients are denoted by fi mn = fi m (ki n )^ 
and K £n is the normalisation constant (see appendices for 
more details). 

In various applications, though, the continuous field / 
cannot be directly observed. This is notably the case in 
cosmology where galaxy surveys give indirect information 



leading to the final Fourier-Bessel coefficients fim(k) 

hm{k) = ^Jdr r 2 f em (r)kj e (kr). (11) 

This formulation hence extends the notion of 2D spherical 
harmonics to three-dimensional fields. 

It is also possible to conceive the reverse approach, i.e. 
to perform the spherical Bessel transform first and subse- 
quently expand the resulting coefficients into spherical har- 
monics. Formally, the ^-th order spherical Bessel transform 
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of / (similar to its Hankel transform) is 
^2 



dk fi(k,0,</))kjt(kr), 



for which the inversion formula gives 



fi(k,o,4>) 



dr r 2 f(r,9,(fi)kji(kr) 



(12) 



(13) 



The result is then expanded into spherical harmonics but 
with an unusual formulation since 0, 0) and Yg m (#, <p) 
(as well as the basis functions ji(kr) and Yg m (0,0)) have 
now the £ parameter in common: 



(14) 



Again, using the inversion formula, we obtain the Fourier- 
Bessel coefficients 



fim(k) 



I <m f £ (k,o,<l>)Y? m (o,< 



(15) 



Due to the closed domains of shells Q(r) and thus the rel- 
ative independence of angular and radial dimensions, the 
raw (equations [I] and |2|, the forward (denoted by SHB 
for spherical- Harmonic^B ess el, equations [8] to [TT| ) and the 
reverse (denote d by SBH for spherical- B ess el^Harmonic, 
equations [l2] to 15 ) methods are equivalent formulations 
of the spherical Fourier-Bessel decomposition of any three- 
dimensional field /(r, #,</>). This is summarised in the fol- 
lowing schematic description of each method: 



RAW : 


m - 


SHB : 


m - 


SBH : 


m - 



three-dimensional integration 



SHT 



fim{r) 



SBT 



> fam(k) 
> flm(k) (16) 



SBT 



SHT 



> fem(k). 



Note that this section is related to the ideal case R = 
oo, but all equations can be straightforwardly rewritten for 
a finite R by replacing k by fc/ n , bounding each integral 
and adapting normalisation. The formulas arising from this 
adaptation are used in the next sections. 

2.3. Estimating Fourier-Bessel coefficients from a real survey 



Although the three approaches described in |2.2| are theoret- 
ically equivalent, their estimates and numerical implemen- 
tations take different forms. 



2.3.1. Forward approach (SHB) 

Estimating the /g mn coefficients using the forward method 
naturally requires the radial dimension to be discretised. 
Indeed, the first step is to compute the spherical har- 
monic transform on a set of shells located at radial values 
n, . . . ,r Nlayers . In each layer, the coefficients fimin) are 
estimated. Although it is possible to perform a raw esti- 
mate for the later harmonics transform, it is often advis- 
able to use a robust 2D discretisation scheme (of N pix (i) 
pixels for the i-th shell) and to take advantage of the re- 
lated high-performance algorithms. Angular space is hence 
discretised into nodes (ri,9 p: (fi p ) = (ri,*y q ) and the field is 



approximated on each node, giving /(r$,7 p ). The spherical 
harmonic decomposition in the i-th shell becomes 



Npix{i) 

hm(ri)= /( r ^7p)^/m(7p), 

p=l 



(17) 



and the final coefficients are obtained by performing the 
following spherical Bessel decomposition: 



hmn = fe m (ri)ke n je(ke n ri). 



(18) 



1=1 



With this method, radial and angular spaces are discretised 
and both transforms are approximated. 

2.3.2. Reverse approach (SBH) 

For the reverse approach, a 2D scheme on the sphere was 
required as well. As previously, this scheme defines a set 
of N P i X zones (pixels) related to angular nodes 7 . If G q 
denotes the points of the survey located in the solid angle 
corresponding to the q-th zone of the scheme, we perform 
the spherical Bessel Transform (raw estimate) in each zone 



finh q) = fl(hn,7q) = fe ^(*^p), 

peG q 



(19) 



and each of these intermediate maps is decomposed into 
spherical harmonic (spherical Harmonics Transform) which 
gives the Fourier-Bessel coefficients 



Npix 
(7=1 



(20) 



With the reverse method, one can avoid to discretise ra- 
dial space. Moreover, this one-shell pixelisation of the sky 
(thus based on physical solid angles) allows for a natural 
treatment of radial distortions (redshift, relativist ic) and 
masking effects. Using multiple resolutions at different ra- 
dial values, as would be possible with the forward method, 
is much more questionable. The SHB method also proves 
to be a powerful tool for weighting the data prior to esti- 
mating the power spectrum. For instance, in Tadros et al. 
( |1999| ) used a fiducial power spectrum to derive an optimal 
weighting operation. This operation is quite complex when 
using the raw Fourier-Bessel approach, whereas the SHB 
formulation naturally handles the dependence on k of the 
weighting function. 

The three methods to estimate the spherical Fourier- 
Bessel decomposition can therefore also be expressed for a 
discrete 3D survey, summarised schematically below: 



^ ATTT r ~~i Raw sum, best estimate of FB coefficients ~ 

RAW : {r p } > f £rnn 

r . Approx SHT 7 / \ Approx SBT ~ 

SHB : {r p \ > fimyn) > fe m n 

OT ^ TT r 1 Exact SBT 7 , \ Approx SHT ~ 

SBH : {r p } > hn\l p ) > ftmn- 

Note that in practice, the range of (Z,ra, n) is finite, 
which introduces an additional approximation. Here, £ and 
n are restricted to [0, imax] and [1, n max ] respectively. Given 
£, m goes from — i to £. 
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3. Method comparison 

3.1. Complexity, accuracy and discretisation grids 

For a survey that probes a field by N discrete points, the 
raw method is the natural estimate of the Fourier-Bessel 
coefficients. However, since each point contributes to the 
calculation of every coefficient fi mn (V /,m,n), computa- 
tion time is proportional to TV • n max (£ ma:c -f l) 2 /2, which 
can be highly problematic for large surveys. 

In the forward method, the repeated spherical har- 
monic transforms take advantage of tesselation schemes 
and high-perf ormance algori t hms s uch as tho se provided 



by HEALPix |G6rski et al.| (|2QQ5|)1, IGLOO |Crittenden| 



fc Turok| ( [1998| )j or GLESP |Doroshkevich et al.| Q2008)j. 

"dere * 



Roughly speaking, the number of nodes to be considered is 
reduced from N to N pixi and the use of fast spherical har- 
monic transforms on these schemes significantly decreases 
computation time. 

However, this approach requires the three-dimensional 
space to be divided into shells ^( r i)- Both radial and an- 
gular dimensions are discretised, and the survey is approxi- 
mated on an actual 3D grid. In practice, this approximation 
deteriorates the accuracy of the estimated Fourier-Bessel 
coefficients. Furthermore, designing a meaningful radial dis- 
cretisation is a difficult task. For equal-area pixelisations, 
the area of each pixel on the z-th shell is Airrf /N pix (i). 
With HEALPix, the riside angular parameter may only be 
increased by a factor 2, which changes the number of pixels 
by a factor of 4 (since N pix (i) — 12n s ^ e (i) 2 ). This means 
that pixel areas cannot be stabilised for subsequent shells 
as r increases. Consequently, it is not possible to adapt a 
resolution to obtain a 3D scheme with equal-volume vox- 
els. Extending 2D schemes to 3D is difficult and may even 
require a novel formalism for an equal- voxel 3D grid. 

In the reverse approach, though, the use of angular 2D 
schemes is possible, but radial space does not need to be 
discretised. The points of the survey are grouped according 
to angular zones instead of being approximated on a 3D 
grid. An estimate of the spherical Bessel transform is com- 
puted in every solid angle, and the result is then expanded 
in spherical harmonics on the 2D spherical grid. In the fi- 
nal account, this method naturally leads to more accurate 
coefficients than the forward method and also takes advan- 
tage of high-performance 2D schemes. For these reasons, 
we focus on the reverse approach and its implementation, 
using HEALPix for the angular transform. 

Finally, for both forward and reverse methods the spher- 
ical harmonics discretised basis (coefficients Yimicivi) ma y 
be fully pre-computed and stored in external files. This is a 
particularly useful feature (incompatible with the original 
formulation of spherical Fourier-Bessel), which significantly 
eases and speeds up the use of these methods. 



3.2. Speed and accuracy of the Reverse Method 

To test the accuracy and speed of the reverse method com- 
pared to the raw method, we considered the high- resolution 
full- sky Horizon simulation ( Teyssier et al.|[2QQ9 ). Horizon 
is a N-body simulation covering a 2/z _i Gpc periodic box 
using 70 billion dark matter p articles using a WMAP3 cos- 
mology (Spergel et al. 2QQ7[ ). A derived halo catalogue is 
available, which we used to calculate fimn and Cg{kgn) val- 
ues using both methods (raw and reverse). Since we are 



interested only in comparing the speed of each method, we 
simply considered each halo to have equal weight. 

We performed the raw and the reverse estimates on 
three 'surveys' of N = 4.2 x 10 5 ,3.1 x 10 6 and 1 x 10 7 
halos, which correspond to three different depths (z max = 
0.1, 0.2 and 0.3 respectively) in the Horizon simulation. The 
HEALPix angular parameter is given by n s ^ e - 

The results of the accuracy and speed tests are given 
in Table 1. The third (fourth) column gives the percent- 
age f coefficients for which the relative accuracy e(/^ mn ) 
(e(CV n )) is lower than 0.3% for given values of n s ^ e and 
N. We considered the intervals (i,n) G ([0,20], [1,20]) and 
(Z, n) G ([21, 50], [21, 50]) separately, since the estimation of 
higher coefficients depends more on the value of n s id e - We 
also compared computation times of the two methods by 
observing the ratio T = t reverse /t raw . Given a survey and a 
method, computation time denotes the CPU time required 
to compute the fc/ n 's (from the Bessel functions) and the 
final coefficients ft mn without using pre-computed quan- 
tities. Note that we performed this analysis by distribut- 
ing the calculations on five machines and simply adding 
the individual contributions to computation time since our 
method is linear with survey size. With the reverse method, 
though, the roots of the Bessel functions as well as the 
spherical harmonics may be pre-computed and stored in 
external files, which decreases computation time and com- 
plexity when working with 3DEX. 



N 


reside 


e r {fimn) < 0.3% 


e r (C en ) < 0.3% 


T 






[0,20] / [21,50] 


[0,20] / [21,50] 




4.2e5 


512 


87% / 42% 


99% / 96% 


8 




1024 


95% / 65% 


99% / 98% 


4 




2048 


99% / 84% 


99% / 99% 


2 


3.1e6 


512 


92% / 50% 


99% / 95% 


10 




1024 


98% / 74% 


100% / 100% 


5 




2048 


99% / 90% 


100% / 100% 


2 


9.7e6 


512 


92% / 50% 


100% / 97% 


12 




1024 


97% / 74% 


100% / 100% 


6 




2048 


99% / 90% 


100% / 100% 


3 



Table 1. Estimation of Fourier-Bessel coefficients: compar- 
ison of the new method, the reverse formulation (equations 
19 and 20 using HEALPix discretisation) with the original, 
raw formulation (equation [7]) . The third (fourth) column 
gives the percentage f coefficients for which the relative ac- 
curacy c(fe mn ) {e{Ci n )) is lower than 0.3% for given values 
of rtside and N. T is the ratio of elapsed times of the two 
methods. 



The reverse method is about an order of magnitude 
faster than the raw method, but this depends on the choice 
of rtside- For n si d e = 1024 almost all fimn coefficients in 
the range [0, 20] (for £ and n) have relative error below 1%, 
and 90% have it below 0.3%, whereas over 99% of C(£, k n ) 
coefficients are accurate to < 0.3%. In the range [20,50], 
the accuracy is somewhat degraded due to the extension 
of the HEALPix formalism to 3D surveys. Indeed, for data 
distributed on the sphere, 3D space is very sparse even for 
large surveys. Increasing n s ide to 1024 or 2048 strongly im- 
proves the accuracy for higher orders £. Note that compar- 
isons for £ > 50 are limited by the amount of time the raw 
method takes. 
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1 




\0 



10 
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30 



40 
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60 
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Fig. 1. Speed results (raw and reverse methods) for in- 
creasing summation limits £ max = n maxi for a survey 
of N = 9.7 x 10 6 halos. Dashed lines are the fitted com- 
plexity curves. The reverse formulation is suitable for 
pre-calculations and external storage of the spherical 
harmonics, which was not performed here but enables 
for additional speed improvements. 



Figures 1 and 2 show the time taken for calculations 
as a function of £ m ax and n max (Figure 1), and as a func- 
tion of number of halos (Figure 2). The boxes correspond 
to the raw method, the circles and diamonds to the reverse 
method with n s id e = 512, 1024 respectively. The dashed 
line corresponds to the general rule that the raw method 
scales as Nnm ax (£ m ax + I) 2 /2, whereas the points are all 
estimated from calculations. With £ max = n max = 100 and 
TV = 9.7e6, the raw decomposition took a few days of cal- 
culations, whereas the reverse method only took 12 hours. 
In our formalism, k max is determined by the choice of R 
for the boundary condition and by the band-limit n max 
for spherical Bessel coefficients. For each multipole £ we 
have k max = k inrnax = qin max /R where qe nmax is the n max - 
the root of the £-th spherical Bessel function. Because R 
is usually imposed by the problem or the data, one must 
increase n max to probe smaller radial scales. In fact, a rea- 
sonable approximation (or even a simple plot) shows that 
Qin ~ {£ + 3n). This observation enabled to predict which 
radial scales are probed and how computation time scales 
with kmax, given that we provide the complexity for £ m ax 
and n max . 

One of the main advantages of the reverse method is 
that it is naturally suited to parallel computing because 
it uses HEALPix fast spherical Harmonics Transform rou- 
tines. All previous tests were performed on a recent com- 
puter (i7 processor, 8Go RAM) and take advantage of 
OpenMP (with four threads). More advanced computing 
means (larger RAM and more processors) significantly de- 
crease calculation time. For example, £ m ax = rimax = 128 
with n S ide — 2048 took about an hour with 128 cores 
and 512Go RAM, whereas computation time for the raw 
method was estimated to several days on the same machine. 
Note that the raw method is also suited to parallelisation: 
galaxies may be treated separately by different threads. In 
all experiments, we took advantage of this property and 



Fig. 2. Speed results (raw and reverse methods) for 
increasing survey size, for £ max = n max = 30. Dashed 
lines are the fitted complexity curves. 



performed both raw and reverse decompositions with four 
threads to perform relevant comparisons between the two. 

In terms of the power spectrum, figures 3 and 4 show 
the relative error between the raw and the reverse methods 
both in mode-mode space (£ — n) and in mode-scale space 
{£ — kin). For this comparison we decomposed a survey of 
N = 4.2 x 10 5 halos with z max = 0.1. Figures in mode- mode 
space naturally differ according to the choice of the bound- 
ary R because the latter determines the discrete radial scale 
spectrum {/^ n }, and hence mode n computed with two dif- 
ferent i?'s corresponds to different /c-scales. When compar- 
ing the results from R = 1000 and R = 2000 in mode-scale 
space, we observe that the boundary condition fixes the ex- 
plored scales. The left column is thus complementary to the 
right column to explore higher values of k. Although figure 
[3] gives information about the final coefficients, figure 4 is 
hence more appropriate to see which scales are probed and 
with what accuracy. 

In view of the £ — k^ n space, we see that no fluctuations 
are observed along the k axis up tofe = 0.03hMpc- 1 .In this 
range, fluctuations occur in £ space, which are accurately 
probed with n s ^ e = 512 until £ = 25 but naturally require 
a more precise scheme for i > 25, k > 0.03hMpc _1 (smaller 
scales in physical space). In conclusion, parameter n s id e (as 
well as R) must be chosen depending on the scales one 
wishes to probe. Figures 3 and 4 provide accuracy results 
that are complementary to Table 1. 



4. The 3DEX library 

The 3DEX library requires the HEALPix package (v2.12 
or later) and the CFITSIO library. 3DEX may either by 
installed with an HEALPix-like procedure (configure and 
make commands) or using CMake. The Fortran modules, 
the 3DEX dynamic library and the related executables will 
be created in the relevant directories (see README file for 
more information). 

In addition to the numerical procedures required to 
compute Fourier-Bessel coefficients, various other routines 
are provided in the library, such as those converting redshift 
to comoving distance, computing spherical Bessel functions 
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Fig. 3. Relative error on the power spectrum in mode-mode space C(l,n). We compare the original 
formulation of spherical Fourier-Bessel decomposition with the reverse formulation, testing n s ^ e = 
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range [-0.3%, +0.3%]. 
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and their zeros, reading and writing 3D structures (fimn 
and C/ n ), or reconstructing radial maps from Fourier-Bessel 
coefficients. 

Three executable programmes are generated: 

— survey2almn performs the spherical Fourier-Bessel de- 
composition (reverse method) of a discrete survey with 
input parameters l ma x, ^max, t and n s ide- Outputs are 
the fimn coefficients and the power spectrum. 

— survey 2 almn -interactive is very similar to the previous 
programme, but converts redshift values into comoving 
distance before performing the spherical Fourier-Bessel 
decomposition. In particular, the routine takes a .txt 
file as input, taking into account parameters on the cos- 
mology and on the decomposition. 

— almn2rmap extracts the fimn coefficients from a FITS 
file and reconstructs the field (HEALPix map) at a given 
radius. Inputs are the resolution, the radius and sum- 
mation limits l max and n max , which allows one to recon- 
struct several maps at different scales and resolutions. 

The corresponding calls are given by the examples below. 

> survey2almn survey_thetaphir.dat almn.fits 
cln.fits 20 20 256 2000.0, 

where survey_thetaphir.dat is a survey with columns 
representing </>, r, and the keywords correspond to 
values of i,n,n S ide,R- The output is both the coefficient 
values (almn.fis) and the Fourier-Bessel spectrum (cln.fits). 

> survey 2almn_inter active parameters.txt, 

where parameters.txt is an external file containing 
input parameters for the survey and the cosmology (which 
allows for more flexible use). Finally, for the map recon- 
struction, we can use: 

> almn2rmap almn.fits map. fits 400.0 256 10 
10 2000.0, 

where the keywords correspond to r macc , n s ^ e , A ^, R- 



5. Conclusion 

High-precision cosmology from galaxy and weak lensing 
surveys will require the analysis of 3D data in spherical co- 
ordinates, a situation for which spherical Fourier-Bessel de- 
composition is most suited. Current methods will be inad- 
equate for future planned cosmological surveys, which will 
provide for example galaxy surveys with billions of galaxies, 
compared to millions today. 

We have reviewed the forward or SHB (spherical 
Harmonic-Bessel) formalism of the spherical Fourier-Bessel 
decomposition which first calculates the tangential, then 
the radial decomposition. We also introduced the reverse 
or SBH (spherical Bessel Harmonic) formalism that inverses 
this order. Only the latter approach can take advantage of 
existing fast codes for the calculation of tangential modes. 
(To do the same, the former would require a a new voxeli- 
sation scheme.) 

We presented a public code 3DEX (3D Expansions) for 
the fast calculation of Fourier-Bessel coefficients and spec- 
tra, which uses the HEALPix pixelisation scheme for cal- 
culating the tangential modes. The 3DEX code is based on 



the reverse /SBH formulation of the Fourier-Bessel decom- 
position. 

We tested the 3DEX code on linear and quasi-linear 
scales (£ < 50 and kt n < 0.2hMpc _1 ) using the Horizon 
halo simulation for redshift s z < 0.3. For n s id e — 1024 the 
3DEX method for calculating the power spectrum C(£,k) 
is accurate to 0.3% on these scales. 

For surveys with < 10 million galaxies, computation 
time is reduced by a factor 4-12 depending on the desired 
scales and accuracy. For larger surveys the gain in time will 
be even greater. Finally, the use of the 3DEX code is not 
restricted to cosmological calculations, and can be used in 
any other discipline that requires a spherical Fourier-Bessel 
analysis of 3D data. 
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Appendix A: Normalisation and discrete radial 
spectrum 

The basis functions kji(kr)Yi m (6, <j>) form a set of eigenfunc- 
tions of the Laplacian operator in spherical coordinates. In 
particular, these functions are orthonormalised in the con- 
tinuous case, thanks to the orthogonality relation 



j dtldr r 2 j l {kr)j l ,{k'r)Y lm {e,cj>)Y l 1 m ,{e,4>) 
~ n 7 5 D (k-k')6«5« m „ 



2kk' 



(A.1) 



(Baddour |2010b where S K is Kronecker's delta notation and 
o D Dirac s function. 

A common approach to simplify the problem is to as- 
sume some boundary conditions for the field /. Different 
conditions have been explored in the literature ( |Fisher et al. 
1995| |Heavens fc Taylor] 1995| ) , including potential or gra- 
dient continuity. In this paper, we used a condition that 
derives from the classical formulation of the discrete spher- 
ical Bessel transform: space is assumed to be finite and 
limited to a sphere of radius R. In this case, the spherical 
Bessel functions are not normalised and the boundary ef- 
fect leads to a discrete spectrum {ki n }. The Fourier-Bessel 
coefficients become a set fimn — flm(kln), and the complete 
description of the field in the so-called Fourier-Bessel basis 
QBinney fc Quinn||1991 ) is summarised in equation [5] 

As a consequence, a natural choice for the boundary 
condition i s to impose the field to vanish at r = R ( jAbramo" 
et al.|20iQ ), which constrains the Bessel functions and gen- 
erates the radial spectrum {ki n } such that, for all I and 
n, 

ji(k ln R) = 0. (A.2) 

If qi n denotes the n-th root of ji(z), the closure relation of 
the Bessel basis is 

JO 



alz z 2 ji(qi n z)ji>(qi> n iz) 



(A.3) 
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which gives with k\ n = qi n /R, 
Jo 



dr r 2 ki n ki> n >ji(ki n r)ji>(ki> n >r) 



tf n \Jl+l(<lln)]' 



■Sll'dnn' 



(A.4) 



The discrete spectrum is thus fixed by the zeros of the 
spherical Besse l functions. We ob tain the normalisation co- 
efficients K{ n ( |Fisher et aT|l995| > 



i? 3 



[ki n ji +1 (ki n R)] 2 , 



(A.5) 



which are used for field reconstruction (equation [5| . 

Other approaches are possible to tackle bounaary con- 
ditions in radial sp ace, notably those imposing potential 
continuity at r = R QFisher et al.||1995 ). Then, the discrete 
spectrum k[ n is such that 



j,_i(fc,' n fl) = 0, 

and normalisation constraint becomes 
-l R 3 

"'in = -J-iklnMhnR)] 2 \ 

Appendix B: Angular masks 



(A.6) 
(A.7) 



3DEX takes into account optional angular masks under the 
form of either an equatorial cut or an input all-sky FITS 
map. 

In the first case, supplying 6 cut defines the latitude (in 
degrees) of a straight symmetric cut around the equator. 
Pixels located within that cut (I = cos(O cut )) are ignored. 

In the second case, the supplied mask must be an 
HEALPix map (ring ordering) of N p i x pixels at resolution 
nside (which must be identical to Fourier-Bessel resolution 
parameter) 



(B.l) 



In the forward method, the first step is to apply the 
mask to each discrete shell before performing the spherical 
Harmonics Transform. Hence for the ith shell, field / is 
weighted by the mask at each pixel 



f\ r i,l q ) = w(<y q )f(ri,<y q ). 



(B.2) 



The Fourier-Bessel coefficients are obtained after perform- 
ing SH and SB transforms. 

In the reverse method, the first step is still the spherical 
Bessel Transform, which gives a set of n max HEALPix maps 
f(kin:Jq)- The mask is then applied to each of these maps 



f'(hn,7q) = ™{lq)f{hn,lq)- 



(B.3) 



and the modified spherical Harmonics Transform gives the 
final fi mn coefficients. 
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